#!/usr/bin/perl -w
use strict;
if($#ARGV + 1 < 3){
        die "Insufficient arguments\nUsage: find_var.pl <pileup-file-name> <wgsim-report-file> <var-report-filename>\n";
        }
my $pileup_fh;
my $report_fh;

open $report_fh,$ARGV[1] or die "can not opem $!\n";
#my @report_cont=<$report_fh>;
#close $report_fh;
my $threshold = 0.2;
my %insertions_report;
my %deletions_report;
my %indel_r;
print STDERR "processing $ARGV[1] ...\n";
while(<$report_fh>){

#foreach my $line (@report_cont){
#	chomp($line);
	chomp;
#	if($line =~ /(.*)(\s+)(\d+)(\s+)([\w-]+)(\s+)(\w+|-)(\s+)(.*)/){
	if($_ =~ /(.*)(\s+)(\d+)(\s+)([\w-]+)(\s+)(\w+|-)(\s+)(.*)/){
		if($5 eq '-'){
			$insertions_report{$3}=$7;
			my $indel_len = scalar(split(//,$7));
			$indel_r{$3} = '+'.$indel_len.$7;
		}
		elsif($7 eq '-'){
			$deletions_report{$3}=$5;
                        my $indel_len = scalar(split(//,$5));
                        $indel_r{$3} = '-'.$indel_len.$5;
		}
	}	
}

print STDERR "processing $ARGV[0] ...\n";

open $pileup_fh,$ARGV[0] or die "can not open $!\n";
#my @pileup_cont = <$pileup_fh>;
#close $pileup_fh;
my %indel_p;
my %cov_p;
while(<$pileup_fh>){
#foreach my $line (@pileup_cont){
#        chomp($line);
	chomp;
	my $line = $_;
	if($line =~ /\+|\-/){
		if($line =~ /(.*)(\s+)(\d+)(\s+)(A|T|G|C)(\s+)(\d+)(\s+)(.*)(\s+)(.*)/){
			my $scan = $9;
			my $ref_num = $3;
			my %count_indel;
			while($scan =~ /(\+|\-)(\d)(\w+)(\W*)(.*)/){
				$count_indel{$1.$2.$3}++;
				$scan = $4;	
			}
			my $max = 0;
			my $max_key;
			my $counter = 0;
			for my $key(sort(keys %count_indel)){
				if($max < $count_indel{$key}){
					$max = $count_indel{$key};
					$max_key = $key;
				}
			#	$max = $key;
			#	$counter++;
			}
			$indel_p{$3} = $max_key;	
			$cov_p{$3} = $max/$7;
		}
	}if($line =~ /\*/){
		
	}
}

print STDERR "comparing...\n";

my $count_indels_report = scalar keys %indel_r;
my @match;
my @diff_fp;
for my $key ( sort (keys %indel_p) ) {
	if(defined $indel_r{$key}){
		#my $temp = "Match at ".$key." $indel_r{$key} \n";
		my $temp = $key." $indel_r{$key} Matches\n";
                push @match, $temp;
	}else{
		if($cov_p{$key} > $threshold){
			#my $temp = "False positive at $key. $indel_p{$key} .Covergae: $cov_p{$key}\n";
			my $temp = "$key. $indel_p{$key} False Positive. Covergae: $cov_p{$key}\n";
			push @diff_fp,$temp;
		}
	}
}
my @diff_fn;
for my $key ( sort (keys %indel_r) ) {
        if(!defined $indel_p{$key}){
		#my $temp = "False negative at $key. $indel_r{$key} .Covergae: 0\n";
		my $temp = "$key. $indel_r{$key} False Negative. Covergae: 0\n";
		push @diff_fn,$temp;
	}
}
my $find_var_fh;
my $find_var_fn = $ARGV[2]; 
open $find_var_fh, ">".$find_var_fn or die("Can not open at $!\n");

my $match_count = scalar(@match);
print STDERR "---------\nMATCH: $match_count\n---------\n";
my $count_line = 0;
foreach my $entry(@match){
	$count_line++;
	#print $find_var_fh $count_line.". ".$entry;
	print $find_var_fh $entry;
}

my $fp_count = scalar(@diff_fp);
print STDERR "---------\nFALSE POSITIVES: $fp_count\n---------\n";
$count_line = 0;
foreach my $entry(@diff_fp){
	$count_line++;
	#print  $find_var_fh  $count_line.". ".$entry;
	print $find_var_fh $entry;
}

my $fn_count = scalar(@diff_fn);
print STDERR "---------\nFALSE NEGATIVES: $fn_count\n---------\n";
$count_line = 0;
foreach my $entry(@diff_fn){
	$count_line++;
        #print  $find_var_fh $count_line.". ".$entry;
	print $find_var_fh $entry;
}
close  $find_var_fh;

